function [ m_ccc,m_css,m_cww,m_cee,m_cnn ] = rrb( Nx , Ny , cc , cs, cw)
%Construction of RRB Preconditioner
%   Detailed explanation goes here

m_ccc =cc;
m_css =cs;
m_cww =cw;

m_cnn = zeros(Nx+3,Ny+3);
m_cee = zeros(Nx+3,Ny+3);
for jj= 1:Ny+2
    for ii=1:Nx+2        
        m_cnn(ii,jj) = cs (ii,jj+1);
        m_cee (ii,jj) = cw (ii+1,jj);
    end
end


m_cne = zeros(Nx+3,Ny+3);
m_cnw = zeros(Nx+3,Ny+3);
m_cse = zeros(Nx+3,Ny+3);
m_csw = zeros(Nx+3,Ny+3);

m_Nx1 = Nx;
m_Nx2 = Ny;

st = 2;
tw = 2;
niv = 0;
while tw < m_Nx2+3
    if (niv > 0)        
       % fac = 1.0 + (st-1) * 10.0 / ((m_Nx1 + 3) * (m_Nx2 + 3));
       fac=1.0;
        for jj = st:tw:m_Nx2+3
            for ii = 1:tw:m_Nx1+3
                m_ccc(ii,jj) = m_ccc(ii,jj) * fac +m_cne(ii,jj) + m_cnw(ii,jj) +m_cse(ii,jj) + m_csw(ii,jj);
            end
        end
        
        for jj =1:tw:m_Nx2+3
            for ii = st:tw:m_Nx1+3
                m_ccc(ii,jj) = m_ccc(ii,jj) * fac +m_cne(ii,jj) + m_cnw(ii,jj) +m_cse(ii,jj) + m_csw(ii,jj);
            end
        end
    end
    
    niv=niv+1
    
    for jj = st:tw:m_Nx2+3
        for ii = st:tw:m_Nx1+3
            iw = ii - st;
            ie = ii + st;
            js = jj - st;
            jn = jj + st;
            if (jn <= m_Nx2 + 3)
                m_cne(ii,jj) = m_cne(ii,jj) - m_cnn(ii,jj) * m_cee(ii,jn) / m_ccc(ii,jn);
                m_cnw(ii,jj) = m_cnw(ii,jj) - m_cnn(ii,jj) * m_cww(ii,jn) / m_ccc(ii,jn);
                m_ccc(ii,jj) = m_ccc(ii,jj) - m_cnn(ii,jj) * m_css(ii,jn) / m_ccc(ii,jn);
                m_cnn(ii,jj) = -m_cnn(ii,jj) * m_cnn(ii,jn) / m_ccc(ii,jn);
            end
            if (js > 1)
                
                m_cse(ii,jj) = m_cse(ii,jj) - m_css(ii,jj) * m_cee(ii,js) / m_ccc(ii,js);
                m_csw(ii,jj) = m_csw(ii,jj) - m_css(ii,jj) * m_cww(ii,js) / m_ccc(ii,js);
                m_ccc(ii,jj) = m_ccc(ii,jj) - m_css(ii,jj) * m_cnn(ii,js) / m_ccc(ii,js);
                m_css(ii,jj) = -m_css(ii,jj) * m_css(ii,js) / m_ccc(ii,js);
            end
            if (ie <= m_Nx1 + 3)
                
                m_cne(ii,jj) = m_cne(ii,jj) - m_cee(ii,jj) * m_cnn(ie,jj) / m_ccc(ie,jj);
                m_cse(ii,jj) = m_cse(ii,jj) - m_cee(ii,jj) * m_css(ie,jj) / m_ccc(ie,jj);
                m_ccc(ii,jj) = m_ccc(ii,jj) - m_cee(ii,jj) * m_cww(ie,jj) / m_ccc(ie,jj);
                m_cee(ii,jj) = -m_cee(ii,jj) * m_cee(ie,jj) / m_ccc(ie,jj);
            end
            if (iw > 1)
                
                m_cnw(ii,jj) = m_cnw(ii,jj) - m_cww(ii,jj) * m_cnn(iw,jj) / m_ccc(iw,jj);
                m_csw(ii,jj) = m_csw(ii,jj) - m_cww(ii,jj) * m_css(iw,jj) / m_ccc(iw,jj);
                m_ccc(ii,jj) = m_ccc(ii,jj) - m_cww(ii,jj) * m_cee(iw,jj) / m_ccc(iw,jj);
                m_cww(ii,jj) = -m_cww(ii,jj) * m_cww(iw,jj) / m_ccc(iw,jj);
            end
        end
    end
    
    for jj = 1:tw:m_Nx2+3
        for ii = 1:tw:m_Nx1+3
            iw = ii - st;
            ie = ii + st;
            js = jj - st;
            jn = jj + st;
            if (jn <= m_Nx2 + 3)
                m_cne(ii,jj) = m_cne(ii,jj) - m_cnn(ii,jj) * m_cee(ii,jn) / m_ccc(ii,jn);
                m_cnw(ii,jj) = m_cnw(ii,jj) - m_cnn(ii,jj) * m_cww(ii,jn) / m_ccc(ii,jn);
                m_ccc(ii,jj) = m_ccc(ii,jj) - m_cnn(ii,jj) * m_css(ii,jn) / m_ccc(ii,jn);
                m_cnn(ii,jj) = -m_cnn(ii,jj) * m_cnn(ii,jn) / m_ccc(ii,jn);
            end
            if (js > 1)
                
                m_cse(ii,jj) = m_cse(ii,jj) - m_css(ii,jj) * m_cee(ii,js) / m_ccc(ii,js);
                m_csw(ii,jj) = m_csw(ii,jj) - m_css(ii,jj) * m_cww(ii,js) / m_ccc(ii,js);
                m_ccc(ii,jj) = m_ccc(ii,jj) - m_css(ii,jj) * m_cnn(ii,js) / m_ccc(ii,js);
                m_css(ii,jj) = -m_css(ii,jj) * m_css(ii,js) / m_ccc(ii,js);
            end
            if (ie <= m_Nx1 + 3)
                
                m_cne(ii,jj) = m_cne(ii,jj) - m_cee(ii,jj) * m_cnn(ie,jj) / m_ccc(ie,jj);
                m_cse(ii,jj) = m_cse(ii,jj) - m_cee(ii,jj) * m_css(ie,jj) / m_ccc(ie,jj);
                m_ccc(ii,jj) = m_ccc(ii,jj) - m_cee(ii,jj) * m_cww(ie,jj) / m_ccc(ie,jj);
                m_cee(ii,jj) = -m_cee(ii,jj) * m_cee(ie,jj) / m_ccc(ie,jj);
            end
            if (iw > 1)
                
                m_cnw(ii,jj) = m_cnw(ii,jj) - m_cww(ii,jj) * m_cnn(iw,jj) / m_ccc(iw,jj);
                m_csw(ii,jj) = m_csw(ii,jj) - m_cww(ii,jj) * m_css(iw,jj) / m_ccc(iw,jj);
                m_ccc(ii,jj) = m_ccc(ii,jj) - m_cww(ii,jj) * m_cee(iw,jj) / m_ccc(iw,jj);
                m_cww(ii,jj) = -m_cww(ii,jj) * m_cww(iw,jj) / m_ccc(iw,jj);
            end
        end
    end
    
    for jj = st:tw:m_Nx2+3
        for ii = st:tw:m_Nx1+3
            m_ccc(ii,jj) = m_ccc(ii,jj) +m_cnn(ii,jj) +m_cee(ii,jj) +m_css(ii,jj) +m_cww(ii,jj);
            
            %                 // Om geheugen te besparen gaan we cnw(ii,jj) opslaan in cnn(ii,jj);
            %                 //                                 cne(ii,jj) opslaan in cee(ii,jj);
            %                 //                                 csw(ii,jj) opslaan in cww(ii,jj);
            %                 //                                 cse(ii,jj) opslaan in css(ii,jj);
            m_cnn(ii,jj) = m_cnw(ii,jj);
            m_cee(ii,jj) = m_cne(ii,jj);
            m_cww(ii,jj) = m_csw(ii,jj);
            m_css(ii,jj) = m_cse(ii,jj);
            
        end
    end
    niv=niv+1;
    
    %
    %             // Uitgaande van het stencil
    %         //                          cnn
    %         //                     cnw      cne
    %         //                 cww      ccc      cee
    %         //                     csw      cse
    %         //                          css
    %         // elimineren we cnw en cne met de scheve 5-punts formule. Dit geeft
    %         //                 cnw      cnn      cne
    %         //                 cww      ccc      cee
    %         //                 csw      css      cse
    %         // We zijn nu weer terug de eerste 9-punts formule, maar op een
    %         // grover rooster.
    for jj = 1:tw:m_Nx2+3
        for ii = 1:tw:m_Nx1+3
            
            iw = ii - st;
            ie = ii + st;
            js = jj - st;
            jn = jj + st;
            if (jn <= m_Nx2 + 3 && ie <= m_Nx1 + 3)
                
                m_cnn(ii,jj) = m_cnn(ii,jj) - m_cne(ii,jj) * m_cnw(ie,jn) / m_ccc(ie,jn);
                m_cee(ii,jj) = m_cee(ii,jj) - m_cne(ii,jj) * m_cse(ie,jn) / m_ccc(ie,jn);
                m_ccc(ii,jj) = m_ccc(ii,jj) - m_cne(ii,jj) * m_csw(ie,jn) / m_ccc(ie,jn);
                m_cne(ii,jj) = -m_cne(ii,jj) * m_cne(ie,jn) / m_ccc(ie,jn);
            end
            if (iw > 1 && jn <= m_Nx2 + 3)
                
                m_cnn(ii,jj) = m_cnn(ii,jj) - m_cnw(ii,jj) * m_cne(iw,jn) / m_ccc(iw,jn);
                m_cww(ii,jj) = m_cww(ii,jj) - m_cnw(ii,jj) * m_csw(iw,jn) / m_ccc(iw,jn);
                m_ccc(ii,jj) = m_ccc(ii,jj) - m_cnw(ii,jj) * m_cse(iw,jn) / m_ccc(iw,jn);
                m_cnw(ii,jj) = -m_cnw(ii,jj) * m_cnw(iw,jn) / m_ccc(iw,jn);
            end
            if (js > 1 && iw > 1)
                
                m_css(ii,jj) = m_css(ii,jj) - m_csw(ii,jj) * m_cse(iw,js) / m_ccc(iw,js);
                m_cww(ii,jj) = m_cww(ii,jj) - m_csw(ii,jj) * m_cnw(iw,js) / m_ccc(iw,js);
                m_ccc(ii,jj) = m_ccc(ii,jj) - m_csw(ii,jj) * m_cne(iw,js) / m_ccc(iw,js);
                m_csw(ii,jj) = -m_csw(ii,jj) * m_csw(iw,js) / m_ccc(iw,js);
            end
            if (js > 1 && ie <= m_Nx1 + 3)
                
                m_css(ii,jj) = m_css(ii,jj) - m_cse(ii,jj) * m_csw(ie,js) / m_ccc(ie,js);
                m_cee(ii,jj) = m_cee(ii,jj) - m_cse(ii,jj) * m_cne(ie,js) / m_ccc(ie,js);
                m_ccc(ii,jj) = m_ccc(ii,jj) - m_cse(ii,jj) * m_cnw(ie,js) / m_ccc(ie,js);
                m_cse(ii,jj) = -m_cse(ii,jj) * m_cse(ie,js) / m_ccc(ie,js);
            end
        end
    end
    st = tw+1
    tw = tw*2
end


end

